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ABSTRACT 

We are carrying out a programme of non-linear time-dependent numerical calcula- 
tions to study the evolution of the thermal instability driven by radiation pressure in 
transonic accretion discs around black holes. In our previous studies we first investi- 
gated the original version of the slim-disc model with low viscosity (a = 0.001) for a 
stellar-mass {IQMq) black hole, comparing the behaviour seen with results from local 
stability analysis (which were broadly confirmed) . In some of the unstable models, we 
saw a violently evolving shock- like feature appearing near to the sonic point. Next, we 
retained the original model simplifications but considered a higher value of a (= 0.1) 
and demonstrated the existence of limit-cycle behaviour under suitable circumstances. 
The present paper describes more elaborate calculations with a more physical viscos- 
ity prescription and including a vertically integrated treatment of acceleration in the 
vertical direction. Limit-cycle behaviour is still found for a model with a = 0.1, giving 
a strong motivation to look for its presence in observational data. 
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INTRODUCTION 



The aim of the research programme, of which the present 
paper forms a part, is to study the origin of the incipient 
instabilities operating in accretion discs around black holes 
and to investigate their observational consequences for di- 
verse astronomical phenomena. For doing this, we use a class 
of simple vertically-integrated, non-self-gravitating models 
of transonic accretion discs. This is well justified because, 
despite the ongoing development of increasingly sophisti- 
cated large-scale numerical models, simple disc modelling 
still remains the central link between theory and observa- 
tions. Some of the results of disc phenomenology might be 
truly fundamental, while others might have more limited do- 
mains of applicability (Balbus & Papaloizou, 1999). It would 
be particularly valuable to determine the sensitivity of the 
disc solutions to the form of the viscous stress and here 
we study how the global evolution of the thermal instabil- 
ity driven by radiation pressure depends on the viscosity 
prescription. This investigation follows on naturally from 
the overall strategy of our research programme. The model 
treated in our previous papers (Szuszkiewicz & Miller 1997, 
1998 - hereafter Papers I and II) did not yet include sev- 
eral important non-local effects, and in particular, it did not 
contain a diffusion-type formulation for the viscosity. Those 
results represent a suitable standard reference for making 



comparison with the results from our present calculations in 
which a diffusion-type formulation for viscosity replaces the 
ap prescription used previously. 

The ipr component of the viscous stress tensor, which 
is responsible for the transport of angular momentum in a 
disc, is normally written in the standard form 



: pvr 



dr ' 



(1) 



where p is the density, i/ is a kinematic viscosity coefficient 
and Jl is the angular velocity of matter in the disc. How- 
ever, the relevant viscosity mechanism for accretion discs is 
clearly not molecular viscosity (which is much too small to 
be of interest in these circumstances) but instead is some 
sort of magnetic or turbulent viscosity. For a general turbu- 
lent viscosity. 



I 



turb^turb ' 



(2) 



where v^^^^ and l^.^^^ are the characteristic circulation veloc- 
ity and length-scale for the turbulent cells (see Kato, Fukue 
& Mineshige 1998), and then it is common to write 

u = aHcs, (3) 

where Ca is the local sound speed, H is the half-thickness of 
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the disc and 



a • 



H 



(4) 



is a quantity which must be less than or comparable to unity 
if the turbulence is subsonic with the turbulent cells being 
confined within the disc. This treatment arises from dis- 
cussions of viscosity associated with hydrodynamic turbu- 
lence (von Weizsacker 1948) but a similar formula can be 
appropriate for magnetic viscosity in some circumstances 
(Shakura & Sunyaev 1973; see also Hawley, Balbus & Stone 
2001). 

The form of the viscosity plays a key role for the way 
in which the disc responds to thermal fluctuations. An in- 
finitesimal perturbation giving a small rise in temperature 
would cause the disc to swell up and, according to Eq. (3), 
would have the tendency to make it more viscous giving 
increased heating (although this also depends on how a 
changes in response to the perturbation). At the same time, 
the cooling rate will also become higher as a result of the 
higher temperature. If the increase of cooling is greater than 
the increase in heating, the situation is stable with the con- 
figuration returning back towards its initial state; however, if 
the increase of heating is greater than the increase in cooling, 
the situation is thermally unstable with the temperature con- 
tinuing to increase away from the initial value until the effect 
is eventually saturated in a non-linear regime. Whether this 
kind of thermal instability can occur in practice for real discs 
depends on how fast the viscosity increases as the temper- 
ature is increased and this is something which is still under 
debate. Our view is that it is certainly of interest to study 
the possible consequences of this thermal instability if it 
were to arise and then to see whether evidence for it is re- 
vealed by observations. This is the line being taken for the 
present programme of work and, for simplicity, we follow the 
precedent of using the viscosity formula (3) with constant a 
(which does allow for the instability to occur if the accretion 
rate is high enough). 

Within the framework of the above discussion, there 
is a further simplification which is frequently made. If one 
considers a stationary Newtonian Keplerian disc and makes 
a number of approximations (including assuming vertical 
hydrostatic equilibrium and integrating analytically in the 
vertical direction), it can be shown that expressions (1) and 
(3) lead to 



Tipr — —ap, 



(5) 



where p is the pressure and a has been rescaled so as to 
remove a numerical constant from the formula. This is the 
well-known ap prescription of Shakura & Sunyaev (1973) 
which has been widely used as a general approximate form 
for the viscous stress even under conditions different from 
those for which the formula was derived. 

Abramowicz et al. (1988), calculated a sequence of 
vertically-integrated slim accretion disc models, using the 
ap viscosity prescription but dropping the Keplerian flow 
condition and including the effects of relativistic gravity in 
an approximate way by using the pseudo-Newtonian poten- 
tial of Paczyriski & Wiita (1980). Doing this, they found 
that if M (the accretion rate) is plotted against the sur- 
face density E for a fixed location in the disc, an S-shaped 
curve is obtained. This suggested the possibility of having a 



limit-cycle behaviour associated with a thermal instability 
driven by radiation pressure. The results of non-linear time- 
dependent calculations performed by Honma, Matsumoto & 
Kato (1991) strongly indicated that such a limit-cycle could 
indeed occur and the first complete solution of this type was 
then subsequently obtained by Szuszkiewicz & Miller (1998). 
This finding, if confirmed by less approximate calculations, 
would have important consequences for the interpretation 
of observations. A first step towards checking its robustness 
is to replace the ap prescription by the one given directly 
by Eq. (1) together with expression (3) for v. Doing this 
makes an important change in the mathematical nature of 
the system of equations, introducing a parabolic diffusion- 
type term into an otherwise hyperbolic system. 



Until now, little work has been directed towards study- 
ing transonic accretion discs with this more general viscosity 
prescription because the direct integration of the equations 
is extremely difficult (Hoshi and Shibazaki 1977, Shibazaki 
1978). The numerical results of Chen & Taam, (1993), who 
examined steady state accretion discs, demonstrated that 
the S-shaped form of the M(S) curve is insensitive to the 
form used for the viscous stress. Papaloizou & Szuszkiewicz 
(1994b) constructed models for the inner part of a transonic 
adiabatic accretion disc, assuming constant specific angu- 
lar momentum and including a full treatment of the verti- 
cal structure, and for comparison purposes constructed the 
corresponding one-dimensional viscous-disc models derived 
under vertical-averaging assumptions (see also Papaloizou & 
Szuszkiewicz, 1994a). They also investigated causality con- 
siderations connected with the formulation of viscosity as a 
diffusion process and concluded that if the specific angular 
momentum gradient is very small when the radial velocities 
are significant, the corrections due to finite propagation ef- 
fects become small so that the Shakura & Sunyaev (1973) 
treatment should be reasonably good. Constraints arising 
from causality considerations may be important though as 
far as the uniqueness of the steady-state flow is concerned. 
Both ap and diffusive viscosity prescriptions have been in- 
vestigated by Artemova, Bisnovatyi-Kogan, Igumenshchev 
and Novikov (2000). They found that the stationary solu- 
tions for the two types of prescription are very similar for 
a < 0.1 but begin to differ at larger a. They also showed that 
the solutions with the diffusive shear stress have one singu- 
lar point which is always of the saddle type. A similar result 
had been obtained by Papaloizou & Szuszkiewicz (1994a) 
for transonic accretion discs with a polytropic equation of 
state. The aim of our present study is to investigate the ef- 
fect of using a diffusive form of viscosity in global non- linear 
time-dependent calculations of thermally unstable slim-disc 
models. 



The plan of the paper is as follows. In Section 2, we dis- 
cuss the set of basic equations used for our present calcula- 
tions, highlighting the changes with respect to the equations 
used in Papers I and II. The modifications in the numerical 
treatment are described in Section 3. In Section 4 we present 
the new results obtained with the revised computer code and 
compare them with those obtained previously with the ap 
prescription. Section 5 contains comments and conclusions. 
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2 WHAT IS NEW IN THE SET OF 
EQUATIONS? 

The basic equations used for describing the evolution of the 
thermal instability in an axisymmetric, non-self-gravitating, 
optically-thick disc were presented in full detail in Paper I 
and the treatment which we will be using for parts of the 
flow which are not optically thick was described in Paper II. 
The assumptions and strategy used for the calculations have 
been fully described in these earlier papers. Here we point 
out and discuss two new features which we are now intro- 
ducing into the set of equations: the direct use of formula 
(1) for the viscous stress (giving dependence on dQ./dr) and 
the introduction of a new dynamical equation for the verti- 
cal acceleration (abandoning the previous approximation of 
taking hydrostatic equilibrium to hold in the vertical direc- 
tion). 

In cylindrical polar coordinates (r, ip, z) centred on 
the black hole (having mass M), our basic hydrodynamical 
equations are as follows. The conservation of mass equation 

and the conservation of radial momentum 

DVr _ Idp ^ jl^-il) 

Dt p dr r'^ 



(7) 



are the same as in Papers I and II. We are using the following 
notation: D/Dt is the Lagrangian derivative given by 

Odd 

E = T.{r,t) is the surface density obtained by vertically in- 
tegrating the mass density p, Vr = Dr/Dt, p is the total 
pressure, I — l{r,t) = rv^{r,t) is the specific angular mo- 
mentum, Ij^ is the value of I for Keplerian motion in the 
pseudo-Newtonian potential with v^p — [GMr / {r ~ r q)"^]^'"^ , 
where r^ is the Schwarzschild radius of the black hole, and 
fijf = Vip/r. Note that Vr is negative for an inflow and that 
the p and p appearing in Eq. (7) refer to values in the equa- 
torial plane. 

The form of the azimuthal equation of motion used in 
the previous calculations 

was obtained from 



Dl_ 
Dt 



_i_d_ 

rS dr 



2 / T^rdz 



(10) 



using the ap viscosity prescription given by Eq. (5). Here, 
instead, we use the expression given directly by Eqs. (1) and 
(3): 



: paCsHr 



on 

dr ' 



and Eq. (10) then gives 

(' 






(11) 



(12) 



Note that, as mentioned earlier, the a appearing here is 
rescaled with respect to that used previously with the ap 
prescription. The value of the rescaling factor depends on 
the way in which the vertical structure is treated. In Paper 



I, we described in detail how this was done in our previous 
work and next we will give a modifled discussion of how it 
is done in the calculations of the present paper. 

The vertical equation of motion can be written in the 
form 

Dvz 
Dt 



p dz dz 



(13) 



where "1> is the pseudo-Newtonian potential given by $ = 
—GM/{R — fg), with R = r + z , and Fzz is a viscous 
force which is set to zero throughout the present work. The 
strategy of the slim disc approach is to proceed by deriving 
a vertically-integrated form of this equation, making use of 
the fact that the disc is taken to be thin enough so that the 
gravitational potential $ can be expanded in the vertical 
direction in powers of z/r and only the lowest-order terms 
retained. We then need two further assumptions about vari- 
ations of quantities in the vertical direction in order to be 
able to make the integration. The pressure and density, p 
and p, are taken to be linked together (in the vertical direc- 
tion) by a polytropic relation p oc p^^^'^ , where A'' is the 
polytropic index. Additionally, we need to assume something 
about the vertical velocity and in the previous work (in com- 
mon with the preceding literature on slim discs) this condi- 
tion was provided by assuming vertical hydrostatic equilib- 
rium to hold so that Dv^/Dt = 0. The equation could then 
be vertically integrated to give (after some manipulation 
and setting N = 2 following Paczyiiski & Bisnovatyi-Kogan 
1981): 

E.. (14) 



nlH^ 



This formula applies for strictly Newtonian gravity as well 
as with the pseudo-Newtonian potential which we are using 
here. Inserting expression (14) into the Newtonian argument 
leading to the ap formula, one obtains that 02 (the a in the 
ap formula) is related to the original a, as in equation (1), 
by 

3^/6 



Q2 



■ai 



(15) 



We will use this relationship for calculating models within 
the new approach which are "equivalent" to ones with a 
particular a in the old formulation. 

In Paper I, we commented that neglecting the vertical 
acceleration may sometimes be a rather poor approximation 
in some parts of the disc, as suggested by two-dimensional 
studies (Papaloizou & Szuszkiewicz 1994b), and that inves- 
tigating the effects of including vertical acceleration in a 
consistent way was one of the important issues to be ad- 
dressed subsequently. Contrary to our original strategy of 
introducing additional modifications one by one, the accel- 
eration term has been included here together with the more 
general form for the viscosity. This was done because we 
found that we could not obtain a numerically stable solu- 
tion with the revised viscosity law unless we also included 
the vertical acceleration. We will discuss this further in Sec- 
tion 5. 

The simplest way of dealing with the acceleration term 
is to replace the assumption of hydrostatic equilibrium with 
that of taking z/H to be constant along fluid flow lines. We 
note that while this is not in accordance with the results 
of two-dimensional hydrodynamic simulations (which reveal 
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a complex circulating flow structure) it is consistent with 
the general picture of the slim disc approach and therefore 
seems appropriate here. Introducing this and then making 
the vertical integration as before, one obtains 



DV, 
Dt 



.P 



GM 



(r — rcY 



(?) 



(16) 



where 14 is the vertical component of velocity for a fluid 
element at the top surface of the disc. 

The introduction of the diffusive form of viscosity re- 
quires changes also in the energy balance equation 



DS 



(17) 



Here T is the temperature, S is the entropy per unit mass, 
Qvis is the rate at which heat is generated by viscous friction, 
now given by 



Qvis = pi' [r — I 



dr ) 



(18) 



and Qrad is the rate at which heat is lost or gained by means 
of radiative energy transfer, here given by 



^rad 



F- 



(19) 



with F~ being the radiative flux per unit area away from the 
disc in the vertical direction (we are neglecting heat conduc- 
tion through the disc in the radial direction) . The vertically- 
integrated form of the energy balance equation may bo writ- 
ten in the following way 



DT 
'Dt 



aHr^T{dn/dr) 



TF~ 



p/p 0.67(12 



10.5/3) 0.67pi/(12 - 

(4 - 3/3) T Dp 
(12 - 10.5/3) p Dt 



10.5/3) 



+ 



(20) 



(Note that there is a typographical error in the equivalent 
equation of Paper I - there equation (17). On the right hand 
side there should be a minus sign in front of the first term.) 
If the medium is optically thick we use the expression 



16o-r'* 

KpH 



(21) 



where a is the Stefan-Boltzmann constant and n is the opac- 
ity. The opacity is approximated by the Kramers formula for 
chemical abundances corresponding to those of Population 
I stars 



K = 0.34(l + 6x lO^^pT"^'^) g"Vm^ 



(22) 



The thermodynamic quantities in the equatorial plane are 
taken to obey the equation of state 

p = kpT + pr (23) 

where p^ is the radiation pressure given by Pr = f r**. 

If the medium is not optically thick, wc follow the ap- 
proach of Hure et al. (1994) and write 



F' 



4aT^ 



'^ + V3 + ^ 



(24) 



where r^ and Tp are the Rosseland and Planck mean optical 
depths (equal to K^^pH and KppH). The expression for the 



radiation pressure corresponding to this F is 

For the Rosseland optical depth, we use the expression 

r„ =0.34S(l + 6 X lO^'^pT^^-'^), (26) 

corresponding to the expression for k given by equation (22), 
while for the Planck optical depth we use 

1 



-X% 



4aT^ v^6r; 

where q^^ is the bremsstrahlung cooling rate given by 
q^^ = 1.24 X W'^^Hp'^T^^^ erg cm"^ s"^ 



(27) 



(28) 



Further discussion of this has been given in Paper II. 

Our aim here is to investigate the extent to which the 
inclusion of the more physical viscosity prescription causes 
differences in non-stationary behaviour with respect to that 
of the "standard" models calculated in our previous papers. 



3 WHAT IS NEW IN THE NUMERICAL 
PROCEDURES? 

The equations presented in the previous section have been 
solved numerically using a modified version of the La- 
grangian hydrodynamics code described in detail in Paper 
I, with the accreting matter being divided into a succession 
of comoving radial zones and with the difference scheme fol- 
lowing the standard pattern for one- dimensional Lagrangian 
hydrodynamics. In this section we discuss the modifications 
which have been necessary in order to make a successful cal- 
culation for a model with the new viscosity prescription and 
including the implementation of the dynamical equation for 
motion in the vertical direction. 

The grid structure was modified from that used previ- 
ously in order to give better resolution where this was most 
needed. The gridding was extensively tested to show that it 
was sufficiently fine to be well within the convergent regime 
and, by experiment, we hope to have found a near optimal 
compromise between accuracy and computational expense. 
The inner edge of the grid was set at r « 2.5rg (when- 
ever the innermost comoving zone was entirely inside 2.5rg , 
it was discarded and the grid was remapped following the 
procedure described in Paper I, using piecewise cubic inter- 
polation). The outer boundary was set at r = lO'^r^, far 
enough away so that no perturbation from the inner regions 
could reach as far as this during the time of the calcula- 
tion. Conditions there were essentially unchanging over this 
timescale, making it easy to impose outer boundary condi- 
tions. For the model described in detail in this paper, we 
used a grid of 400 comoving zones organised as follows: for 
the first 59 zones, each one contained a mass 11 per cent 
larger than the one interior to it; for the next 60 zones, the 
increase was 5 per cent; for the rest of the disc, with excep- 
tion of the last 60 zones, the increase was 2.5 per cent and, 
finally, the last 60 zones had an increase of 12 per cent. This 
was found to give satisfactory resolution and accuracy. 

The equations are solved for E, Vr, I, Vz and T, as the 
main dependent variables, with r and H being calculated 
from Dr/Dt = Vr and DH /Dt — Vz respectively. In our 
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previous studies the time-step was adjusted in accordance 
with the Courant condition and with two additional con- 
straints on the fractional variations of p and T in any single 
time-step. Now that we have a parabolic term appearing in 
Eq. (12), it is necessary to add a further time-step condition 
in order to maintain numerical stability: 



At < {Arf/i2u), 



(29) 



where At and Ar are the time-step and space-step respec- 
tively. In the energy equation (20) , which is solved implicitly 
for T, the gradient of the angular velocity appearing in the 
viscous heating term is no longer substituted by the Keple- 
rian value as before. 

When r, S, T and H are known, p can be found directly 
from p = 'E/H while p follows from the equation of state 
(23). Previously, the solution for p, p and H was embed- 
ded in the iteration loop for T. Introducing the dynamical 
equations for 14 and H significantly speeds up the code. 

At every time-step, it is necessary to supply boundary 
conditions for Vr and I at both the inner and outer edges 
of the grid. As before, the boundary conditions at the inner 
edge are set by noting that I is essentially unchanging for 
infalling matter there, so that it can be taken as constant in 
time for the innermost zone, and that the pressure gradient 
term is having negligible effect there, so that it can be set to 
zero when calculating Vr at the inner edge. As noted above, 
the outer edge is located at such a large radius that con- 
ditions there are essentially unchanging over the timescale 
of the calculation and so Vr is kept constant in time there 
while I is calculated using the Keplerian condition at that 
point. 

As initial conditions for the evolutionary calculations, 
we again used a stationary transonic disc model calculated 
with the ap viscosity prescription. The data from this were 
then transferred onto the finite-difference grid used for the 
evolutionary calculations and the a was suitably rescaled 
following Eq. (15) for use in expression (3). The quantities 
transferred were r, Vr, S, T and H; I was then calculated on 
the grid, using the stationary version of equation (7), and 
Vz was calculated from the formula: 



Vz 



dH 

dr 



(30) 



While treating the vertical equation of motion (13) by 
means of the approximate form (16) is a great improvement 
over assuming hydrostatic equilibrium in the vertical direc- 
tion, it is still greatly simplified, particularly in that it in- 
volves taking Dvz/Dt to depend linearly on z/H (which fol- 
lows from z/H being constant along flowlines). This ceases 
to be reasonable when there are rapid dynamic changes in 
the vertical direction such as those seen when the inner parts 
of the disc deflate rapidly in the course of the evolution 
which we will be describing in the next section. In reaching 
a new steady state after the deflation, the vertical accelera- 
tion must inevitably deviate from Dvz/Dt oc z/H. In order 
to keep the solution regular in these circumstances, we intro- 
duced an artificial viscous term Qz oc pVz which augments 
the pressure in the vertical direction when the vertical com- 
pression is very rapid. This, to some extent, mimics the real 
behaviour which would occur in these circumstances and is 
probably the best that can be done within a vertically inte- 
grated approach although, of course, the real answer is to go 



to a full 2D (or 3D) calculation! Related with this problem is 
another one concerned with treating the stage in the imme- 
diate aftermath of the deflation when a sharp density spike 
appears. The rapid compression in the radial direction leads 
to a bunching of the Lagrangian zones causing some of them 
to become so narrow that accuracy is lost. This was com- 
batted by adjusting the regridding routine to avoid having 
ultra-narrow zones and simultaneously increasing the artifi- 
cial diffusion applied in the radial direction. The increased 
value of this artificial diffusion coefficient is needed only for 
a very short time and, while its effect is in the direction of 
causing a broadening of the spike, it is not the main cause 
of the broadening seen which takes place over a much longer 
time than that for which the diffusion is applied. 

As noted in Paper I, the process of transferring the ini- 
tial data onto the finite-difference grid, together with the 
action of numerical noise which inevitably enters during the 
calculation, is sufficient for introducing suitable generalised 
perturbations in the model to trigger growth of unstable 
modes which may be present. 

4 LIMIT CYCLE BEHAVIOUR 

In Paper II we presented a complete limit-cycle solution for 
a model with black hole mass M — WMq, initial accretion 
rate M = O.OGMcr and viscosity parameter a = 0.1. {Mcr 
is the critical accretion rate corresponding to the Eddington 
luminosity). Here we investigate an equivalent model using 
the more physical viscosity prescription given directly by 
equations (1) and (3) (for which T^p oc dfl/dr) and including 
the effects of vertical acceleration by means of Eq. (16). The 
initial equilibrium model is identical to that used in Paper II 
which, according to the local stability criterion, is thermally 
unstable in a region extending from 4.5 r^ to 17.5 r^. 

The main conclusion from our new calculations is that 
the overall behaviour is very similar to that seen before (and 
described in Paper II). In particular, we again see a succes- 
sion of limit-cycles. We here present some additional results 
and highlight some particular points of interest. The figures 
presented here should be viewed in conjunction with those 
of Paper II. 

The evolution of the geometrical profile of the disc is 
shown in Fig. 1 and Figs. 2-5 show the behaviour of the 
temperature T, the surface density S, the effective optical 
depth T^,, and the local accretion rate rh at times corre- 
sponding to the first five panels of Fig. 1. (No curves are 
drawn corresponding to the sixth panel in order to avoid 
the figures becoming over-complicated.) In order to check 
for possible long-term variations, we continued the calcula- 
tion for a total of twenty complete cycles; the results shown 
in Fig. 1 are for a representative cycle - the 17th. 

The first panel of Fig. 1 shows the disc just before the 
start of the cycle (which we define to be the moment when 
the instability first begins to exponentiate). As the insta- 
bility sets in, the temperature in the unstable region rises 
rapidly above its stationary value, increasing the contribu- 
tion of the radiation pressure by nearly an order of magni- 
tude. The heated material expands in all directions, pushing 
inner material into the black hole, expanding the disc verti- 
cally and launching a compression wave out through the disc 
which sweeps matter ahead of it leaving a semi-evacuated re- 
gion behind. This is the stage reached by 2 seconds into the 
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Figure 2. The temperature T, measured in degrees kelvin, is 
plotted against r/r^ at successive times during one full cycle. 
The short-dashed curve shows the temperature at the beginning 
of the cycle (t = Os), the dotted curve is for t = 2s, the dashed 
curve is for t = 8 s, the dot-dashed curve is for t = 16 s and the 
solid curve is for t = 18 s. 
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Figure 1. The evolution of the geometrical profile of the disc 
during one full cycle. The sequence of panels shows the situation 
after 0, 2, 8, 16, 18, 22 and 787 s respectively, starting from the 
beginning of the cycle. (Note that the outer boundary of the grid 
is at r = 10'' r^^ , far beyond the outer edge of the frame shown 
here.) 



cycle, corresponding to the second panel in Fig. 1. At the 
start of the cycle, the effective optical depth t_. , , is greater 
than 100 everywhere but it becomes smaller in the inner 
parts of the disc as the material there expands, falling be- 
low 10 for r < 9rg by the present stage. The behaviour of 
this and the other quantities is shown by the corresponding 
curves (marked with dots) in Figs. 2-5. (Note that despite 
their sharp appearance, all of these curves are, in fact, per- 
fectly continuous and well-resolved on the grid as can be seen 
if they are viewed with a larger spatial resolution.) As the 
outgoing wave progresses further, the temperature peak is 
progressively reduced in magnitude although it remains at a 
level well above that of the initial model. The outgoing wave 
is heating the material through which it passes, causing the 
part of the disc internal to it to swell up further as shown 




Figure 3. The surface density E, measured in units of gem , 
is plotted against r/r^^ at the same times as shown in Fig. 2. 

in the third panel of Fig. 1. The region with very low effec- 
tive optical depth ("r^f , < 10) is now extending out to 15 r,,. 
The compression wave is pushing material strongly outwards 
away from the black hole (note, in Fig. 5, the large negative 
m associated with it) while, behind it, matter is falling in- 
wards at a rate which is far above that in the stationary state 
(m — 0.06) and is even well above the critical accretion rate 
m = 1. In this initial phase of the evolution, the propagat- 
ing compression front divides the disc into two parts. The 



© RAS, MNRAS 



Non-linear evolution of thermally unstable slim discs 7 




Figure 4. The effective optical depth t^,, is plotted against r/r^^ 
at the same times as shown in Fig. 2. 
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Figure 5. The local accretion rate rn(r), measured in units of the 
critical accretion rate Mcr, is plotted against r/r^^ at the same 
times as shown in Fig 

(m = —2nrT,Vr). 



2. Negative values of m signify an outflow, 



inner part has a toroidal structure and is characterised by 
high temperature, low density, r^,, < 100, and a nearly- 
constant supercritical accretion rate, while the outer part is 
basically similar to the original unperturbed stationary disc 
(although, as we will see later, it is actually far from being 
a stationary flow for cycles after the first one) . The width of 
the transition zone between these two parts increases as the 
front moves outwards and has reached 20 r^ by 8 seconds 
into the cycle, corresponding to the third panel of Fig. 1. 

In order for the outgoing wave to continue its prop- 
agation, it is necessary that the material behind it should 



remain in a hot state. Once the wavefront has moved beyond 
the linearly-unstable region, it becomes progressively harder 
for the temperature to be maintained in the high state and 
eventually the front starts to weaken, the temperature drops 
in the innermost part of the disc and the material there de- 
flates, collapsing down into the equatorial plane. Panel 4 of 
Fig. 1 shows the situation at the moment when this happens, 
16 s into the cycle, and in the subsequent panels 5 and 6 one 
can see the deflation spreading out through the disc. During 
the initial collapse, the volume density p rises dramatically 
at the place where the collapse occurs, growing by a factor 
of 10^ in less than a second and forming a sharp spike which 
oscillates and then starts to broaden progressively. (This is 
a feature which is easy to follow in animations of the data 
but is hard to show clearly with a static graph.) At the time 
of panel 4, the accretion rate is decreasing towards the black 
hole and the optical depth has increased considerably in the 
inner parts of the disc. 

Following this rapid change, the temperature continues 
to decrease and drops to a low state well below that of the 
initial model. By the time corresponding to panel 5 of Fig. 1 
(18 s after the beginning of the cycle), the compression front 
is seriously weakening and a growing region of the inner part 
of the disc has deflated. As can be seen from Fig. 3, there 
is very little matter inward of the compression front at this 
stage but the process of refilling the inner part of the disc 
is now beginning, as the accretion flow near to the black 
hole is almost turned off. The disc now consists of three 
different parts. The innermost part forms a geometrically 
thin disc-like configuration with the temperature and surface 
density being lower than those of the initial model and with a 
very low accretion rate. The size of this part is progressively 
increasing; at 18 s it extends out to just beyond 30 r^ . The 
middle part is what remains of the toroidal structure. This 
is still at a higher temperature than in the initial model, 
is under-dense, has a lowered optical depth and has a high 
inward accretion rate. The outermost part is still in a similar 
state to that of the initial model but has been compressed 
somewhat by the still outward- moving compression front. 

At the time of frame 6 of Fig. 1 (22 s after the begin- 
ning of the cycle), the first section of the disc (geometrically 
thin and under-dense) extends out to 65 r^ , the compression 
front has almost died and the remaining toroidal structure 
is shrinking away. 

This is the end of the "outburst" part of the cycle 
(which has proceeded on the thermal timescale) . What hap- 
pens next is a slow and uneventful refilling (and reheating) 
of the inner part of the disc which occurs on the much longer 
viscous timescale. Finally (787 s after the beginning of the 
cycle - Fig. 1, 7th panel) the disc has returned to essentially 
the same state as it was in at the beginning of the cycle. 
The thermal instability then appears again and a new cycle 
is initiated. 

We have calculated the time evolution of the model up 
to the beginning of the twenty-first cycle. The cycle period 
(~ 787 s) is slightly longer than the 780 s found in our pre- 
vious study, and it was seen to be gradually increasing with 
time (albeit with extremely small fractional changes) until 
a fully relaxed solution was finally obtained. In connection 
with this, it is important to stress that after the initial time, 
at which data for a stationary solution was specified, the 
configuration never again passed through a stationary state 
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Figure 6. The local accretion rate profile at the beginning of 
successive cycles (note that the vertical scale is very different 
from that of Fig. 5). The convergence process is rather slow; the 
final converged profile shows rh varying significantly with r and 
so represents a genuinely non-stationary solution. 



and was not even close to doing so. Figure 6 shows the pro- 
file of m at the beginning of successive cycles (just before 
the instability starts) with the initial model being shown 
with the dashed line and the subsequent curves being in or- 
der of progressively decreasing amplitude. Note that for the 
fully relaxed case, rh has a minimum of 0.045 (near to the 
black hole) and a maximum of 0.073 (just outside 100 r^) 
and so is far from being constant at 0.06 as was the case for 
the initial model. In view of this, it is clear that a number 
of cycles would be needed before a relaxed solution could 
be obtained and this is the main reason for the progressive 
change seen in these rh curves although there is also some 
evidence of smoothing related to the parabolic prescription 
for the viscosity. 

Integrating the radiated flux per unit area over the disc 
at successive times, we have obtained the bolometric light 
curve shown in Figure 7. The disc luminosity exhibits a 
burst-like time variation with a burst duration of about 20 s 
and a quiescent phase lasting for the remaining 767 s of the 
cycle. The amplitude of the variation is around two orders 
of magnitude: at the maximum, the luminosity is approxi- 
mately 40 times larger than it was immediately before the 
outburst and after the peak it then drops below the pre- 
burst value by a factor of about 4.5. During the following 
quiescent phase, it then gradually increases again until the 
onset of the next outburst. 

The bolometric luminosity is not a directly observed 
quantity, and so even if we expect that most of the energy in 
galactic black hole candidates will be radiated in the X-ray 
band, it is necessary to calculate the spectrum emitted from 
the disc and the light curves in the particular wavebands, in 
order to perform any detailed comparison with observations. 
This will be discussed in subsequent papers. 



Figure 7. Variation of the bolometric luminosity of the disc ex- 
pressed in units of the Eddington luminosity. 



5 DISCUSSION AND CONCLUSIONS 

In the previous calculations presented in Papers I and II, 
which we use here for comparison with the results of the 
present more elaborate study, we investigated the simple 
slim-disc model with all of its associated assumptions and 
approximations. Two different values of the viscosity param- 
eter Q (low - 0.001 and high - 0.1) were considered there. 
In the present work, two important new aspects have been 
introduced: the diffusion-type formulation for the viscosity 
and the dynamical treatment of motion in the vertical di- 
rection. We now examine some features discussed before for 
the previous formulation to see how these are changed with 
the new one. 

In Paper I we observed that the supposedly stable model 
with a = 0.001 and L = 0.01 L^ did not remain unchanging 
but instead developed an instability in the region 3 — 6 r^ 
after 10 — 20s of the evolution (see Figure 11 in Paper I). 
This appeared first at around 4.7 — 4.8 r^ and then grew 
rapidly, causing termination of the calculation. The high-a 
model of Paper II was subject to the same kind of problem. 
We found that the instability could be suppressed by adding 
a low level of numerical diffusion to the radial equation of 
motion and this suggested that additional physical diffusive 
terms, which were not included in the simple version of the 
model, might play a key role in avoiding this type of instabil- 
ity. Following our strategy, we repeated the calculation for 
the thermally-stable model studied in Paper I, using the new 
version of the code and leaving out the numerical diffusion 
term. We found that the oscillations do appear again (see 
Figure 8) but this time they remain bounded at a low level 
and do not cause termination of the evolution. In particu- 
lar, we note that they are not trapped around their point 
of origin but are able to propagate freely outwards. In the 
case of the high-a model studied here, these oscillations in- 
terfere with the developing thermal instability which is, of 
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Figure 8. Numerical noise which appears when the artificial dif- 
fusion term is removed. These velocity curves are for the model 
with a = 0.001 and m = 0.01 at three different times: t = Os 
(dashed line), t = 76s (dotted line) and t = 1881s (solid line). 
Note that the vertical scale is very much expanded with respect 
to Fig. 11 of Paper 1 so that the oscillations can be seen. 

course, undesirable and so having studied their behaviour 
under these new circumstances, we then re-introduced the 
low level of artificial diffusion as before in order to suppress 
them. 

Next we consider the violently unstable feature seen 
in the calculations presented in Paper I for models with 
low viscosity and luminosities between 0.09 L^ and IL^. 
This instability manifested itself as a dramatic shock-like 
feature which appeared near to the sonic point and then 
grew rapidly, destroying the calculation. Similar behaviour 
was seen in equivalent calculations made by Ryoji Mat- 
sumoto (private communication) using a different numerical 
method. The nature of this behaviour is still under study. If 
not avoided or saturated, it would lead to disruption of the 
disc meaning that models like these would not be suitable 
for describing persistent astronomical flow configurations. 
We made a detailed investigation of the effect of introducing 
different forms of artificial viscosity or diffusion and found 
that none of these could act to limit the instability unless 
they were increased to levels which did not seem reason- 
able. We emphasise that our numerical code has no diffi- 
culty in dealing with ordinary shock structures and, indeed, 
in Paper II we showed a disc evolution with strong shocks 
present in the supersonic part of the flow. In fact, despite 
first impressions, the dynamically growing feature turns out 
not actually to be a shock but, instead, is a transient feature 
across which the Rankine-Hugoniot conditions are not satis- 
fied. It is therefore not particularly surprising that methods 
designed for handling shocks were not successful in deal- 
ing with this. By adding a large enough amount of viscous 
dissipation, the disruptive nature of the feature can indeed 
be removed but our feeling was that introducing so much 
diffusion without physical motivation was unwarranted (as 
we said in Section 4 of Paper II) . At that stage, we wanted 



to wait for the introduction of a (dQ/dr) viscosity to see 
the effect of that before proceeding further. Perhaps a good 
description of the feature is that it is like a wave washing 
up against a sea-wall with the density gradient near to the 
sonic point playing the role of the wall. When sufficiently 
high artificial diffusion was added, the wave was seen to fall 
back from the wall and then to wash up against it again 
repeatedly but without disrupting the rest of the solution. 
The introduction of the additional features in our current 
code has not turned out to be sufficient, in itself, to avoid 
the instability destroying the solution unless similarly high 
levels of artificial diffusion continue to be added. 

Our original plan was to perform a systematic study of 
disc evolution with additional effects being added one at a 
time but here we have gone against this by introducing, at 
the same time, the diffusive (dQ/dr) viscosity prescription 
and the dynamical treatment of motion in the vertical di- 
rection. The reason for this was because of a computational 
difficulty which we encountered with the new viscosity pre- 
scription at the stage when the vertically-expanded part of 
the disc started to deflate. All of our attempts to calculate 
through this stage with the dfl/dr viscosity were unsuc- 
cessful until we decided to introduce also the vertical accel- 
eration, at which point the problem disappeared. It seems 
that the approximate nature of the ap prescription somehow 
compensated the similarly doubtful use of the condition of 
vertical hydrostatic equilibrium. 

In our calculations, we have reported two main types of 
time-varying behaviour for models undergoing the thermal 
instability: for our higher value of a (= 0.1), we have seen 
limit-cycles, while for low a (= 0.001), we saw a seemingly 
catastrophic instability. Takeuchi & Mineshige (1998) have 
presented extremely interesting results from an ap calcula- 
tion for very high a {— 1) in which they see an evolution 
which starts in the same way as our model which makes 
limit cycles but, after the general deflation, a small inflated, 
optically thin region remains at the inner edge of the disc (in- 
ward of r = 5 Tq ) which seems to be a persistent ADAF-like 
feature. This blocks the possibility of starting a new outburst 
cycle. The authors give a persuasive argument for why this 
should occur, indicating that there should be a difference in 
behaviour for models with a > 0.3, for which the expanded 
states are fully optically thin (our model with a = 0.1 has 
considerably reduced optical depth in the expanded region 
but only has r_. , , < 1 briefly and then only for a very small 
region at the inner edge of the disc). When we have, in 
the past, tried to make calculations for q = 1, we experi- 
enced difficulty in keeping the code numerically stable (and 
there is, perhaps, evidence of some similar problem in the 
graphs shown by Takeuchi & Mineshige, although we stress 
that we flnd their results believable, particularly in view of 
the supporting arguments which they present) . Clearly, it is 
now of interest for us to return to those calculations and see 
whether we can confirm the results of Takeuchi & Mineshige. 
In particular, it will be interesting to see whether the per- 
sistent ADAF feature would survive changing the viscosity 
prescription to the diffusive form. 

In conclusion: the main result of our present paper con- 
cerns the effect of replacing the ap viscosity prescription 
by the more physical diffusive form when calculating the 
global evolution of the thermal instability, driven by radi- 
ation pressure, for a vertically-integrated model of a non- 
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self-gravitating transonic accretion disc in the high-a case 
(q = 0.1). We find tfiat the evolution remains cyclic in char- 
acter and is very little changed from that seen before. This 
provides some strengthened motivation for saying that the 
type of disc behaviour seen in our calculations might truly 
have some fundamental importance. It is then natural to 
proceed with studying the possibility that its effects might 
be present among the range of different variability patterns 
observed for accreting compact objects. The transonic discs 
described here with accretion rates larger than a certain 
value (which depends on the mass of the central black hole, 
M, and the viscosity parameter, a) are not stationary but 
show a cyclic behaviour which can be characterised by the 
period of the cycles, the durations of the bursting phase and 
of the quiescence, and the amplitude and shape of the burst. 
For different input parameters: M, a, rh, we should expect 
to see different characteristics. A straightforward compari- 
son with XTE observations, for example, can lead to identi- 
fication of sources where the thermal limit-cycle behaviour 
may be responsible for observed variabilities. So far, there 
is only one source for which this mechanism has been pro- 
posed: GRS 1915-f 105 (Belloni et al. 1997) and the situation 
for that is still under investigation. Good quality data are 
now available for making a detailed confrontation with clear 
model predictions. 



This paper has been produced using the Royal Astronomical 
Society/Blackwell Science TpjjX macros. 
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